import numeric from "numeric"
import _ from "lodash"

/**
 * CPSO算法的javascript实现
 * @author 老朱, 小熊
 * @date 2020-07-04
 */

function CPSO(interNum, dimension, ub, lb, fitness_function, particle_num, maxStep, fit_func_caller = null) {
  // 每个群的粒子个数
  this.particle_num = particle_num;
  // 单个粒子的维度
  this.dimension = dimension;
  // 种群数
  this.group_num = interNum;

  // 单个粒子的每个维度的上限
  this.ub = ub;
  // 单个粒子的每个维度的下限
  this.lb = lb;
  // 单个粒子的适应度函数
  this.fitness_func = fitness_function;
  // 适应度函数的调用者
  this.fit_func_caller = fit_func_caller;
  // 粒子群算法中的粒子列表
  this.particle_list = [];

  // 种群列表
  this.group_list = [];

  // 全局粒子的最优位置
  this.gbest = null;
  // 全局粒子最优的适应度值
  this.gbestValue = Infinity;
  // 当前步数
  this.currentStep = 0;

  this.maxStep = maxStep;
  // 上下文环境变量
  this.b = [];


  /**
   * 粒子群整体初始化
   */
  this.init = function () {
    console.log(`种群数量为:${this.group_num}`)
    // 更新种群
    for (let j = 0; j < this.group_num; j++) {
      // 每个种群增加了id, 表示是对第几个干扰物进行优化
      let group = new Group(j, this.particle_num, this.ub, this.lb, this.evaluate, this.dimension);
      group.init();
      this.group_list.push(group);
      this.b.push(group.groupGbest);
    }
  };

  this.calcInertiaWeight = function() {
    // 惯性权重的两个值
    const w0 = 0.9;
    const w1 = 0.4;
    // 根据当前迭代数和最大迭代数得到当前惯性权重
    //const currentW = w0 - ((w0 - w1) * this.currentStep / this.maxStep);

    //非线性变化
    const currentW = w0 - ((w0 - w1) * (this.currentStep / this.maxStep)**(1/Math.PI/Math.PI));

    //const currentW = (2/this.currentStep)**0.3;

    return currentW;
  }

  this.calcC1 = function() {
    const cStart = 2.5;//1
    const cEnd = 0.5;//3
    const currentC1 = cStart + ((cEnd - cStart) * this.currentStep / this.maxStep);
    return currentC1;
  }
  this.calcC2 = function() {
    const cStart = 0.5;//3
    const cEnd = 2.5;//1
    const currentC2 = cStart + ((cEnd - cStart) * this.currentStep / this.maxStep);
    return currentC2;
  }

  this.step_run = function (steps) {
    // fixed: 通过groupList来判断是否第一次运行, 是的话就初始化粒子群
    if (this.group_list.length === 0) {
      console.log("初始化CPSO")
      this.init();
    }
    // 按照指定的步数更新PSO算法
    for (let i = 0; i < steps; i++) {
      this.currentStep++;  // 步数加1

      // 计算惯性权重
      const w = this.calcInertiaWeight();
      const c1 = this.calcC1();
      const c2 = this.calcC2();

      for (let j = 0; j < this.group_list.length; j++) {
        const group = this.group_list[j]  // 获得种群
        // fixed: 更新的时候j值没传
        group.update(w,c1,c2);//w=0.8,c1=c2=2.0或者1.494
      }
      this.contextUpdate()  // 更新环境变量
    }

    //打印
    console.log(`迭代${this.currentStep}次, 当前全局最优值为${this.gbestValue}`);
    return {
      currentStep: this.currentStep,
      gbest: this.gbest,
      gbestValue: this.gbestValue
    };
  };

  // 更新环境变量
  this.contextUpdate = function () {
    this.b = []  // 更新环境变量
    for (let i = 0; i < this.group_list.length; i++) {
      const group = this.group_list[i]  // 获得子群
      this.b.push(group.groupGbest)
    }
  }



  /**
   * 评估函数, 传入id和要替换的向量值即可求出适应度函数值
   * 使用箭头函数保证this绑定
   */
  this.evaluate = (id, replacePos) => {
    let newContextVec = _.cloneDeep(this.b) // 复制一份环境向量
    if (id && replacePos) {
      newContextVec[id] = replacePos
    }

    // fixed: 把newContextVec拉平, 因为newContextVec是一个一个二维向量的集合
    const flatList = []
    for (let i = 0; i < newContextVec.length; i++) {
      flatList.push(newContextVec[i][0])
      flatList.push(newContextVec[i][1])
    }
    const fitnessValue = this.fitness_func(flatList)
    if (fitnessValue < this.gbestValue) {
      // 如果小于gbest, 则替换更新
      this.gbestValue = fitnessValue
      this.gbest = flatList
    }
    return fitnessValue;
  }


}


function Group(id, particle_num, ub, lb, evaluate, dimension) {
  this.id = id  // 对应维度id
  this.particle_num = particle_num;  // 粒子数目
  this.particle_list = [];  // 存放粒子的数组
  // 单个粒子的每个维度的上限
  this.ub = ub;
  // 单个粒子的每个维度的下限
  this.lb = lb;
  // 替换id和值的评估函数
  this.evaluate = evaluate;
  // 单个粒子的最优位置
  this.groupPbest = null;
  // 单个粒子的最优适应度值，初始值为正无穷大
  this.groupPbestValue = Infinity;
  // 群中粒子的最优位置
  this.groupGbest = null;
  // 群中粒子最优的适应度值
  this.groupGbestValue = Infinity;
  this.dimension = dimension;



  this.init = function () {
    console.log(`种群内粒子数量为:${this.particle_num}`);

    for (let i = 0; i < this.particle_num; i++) {
      let particle = new Particle(this.id, this.dimension, this.ub, this.lb, this.evaluate);
      particle.init();
      this.particle_list.push(particle);

    }
    // fixed: 位置属性应该是position, 而不是pos, 对于群来说,没有pbest, 只有gbest
    this.groupGbest = this.particle_list[0].position;
  };

  // group的更新函数
  this.update = function (w,c1,c2) {
    for (let i = 0; i < this.particle_num; i++) {
      const particle = this.particle_list[i]
      // 获得粒子的适应度值
      const fitnessValue = particle.getFitnessValue()
      // 如果适应度值小于群体中的最优值,则进行替换
      if (fitnessValue < this.groupGbestValue) {
        this.groupGbestValue = fitnessValue
        this.groupGbest = particle.position
      }
    }
    for (let i = 0; i < this.particle_num; i++) {
      this.particle_list[i].update(w,c1,c2,this.groupGbest);//0.4,0.8,0.7298
    }
  }
}

function Particle(id, dimension, ub, lb, evaluate) {
  this.id = id  // 粒子id, 表示优化的第几个
  // 单个粒子的维度
  this.dimension = dimension;
  // 单个粒子的每个维度的上限
  this.ub = ub;
  // 单个粒子的每个维度的下限
  this.lb = lb;
  // 单个粒子的评估函数, 需要传入id和粒子位置
  this.evaluate = evaluate;
  // 当前粒子的位置
  this.position = [];
  // 当前粒子的速度
  this.velocity = [];
  // 单个粒子的最优位置
  this.pbest = null;
  // 单个粒子的最优适应度值，初始值为正无穷大
  this.pbestValue = Infinity;


  // 粒子初始化公式
  this.init = function () {
    // 随机初始化位置
    this.position = this.get_rand_vector();
    // 随机初始化速度,全为0
    this.velocity = numeric.mul(numeric.random([this.dimension]), 0);
    // 当前粒子的pbest就是这个位置
    this.pbest = this.position;
    // 当前粒子的best适应度值就是该位置的评价值
    this.pbestValue = this.evaluate(this.id, this.pbest);

  };

  // 获得在lb和ub范围内的随机向量
  this.get_rand_vector = function () {
    // 先产生一组0-1随机值
    let rand_vec = numeric.random([this.dimension]);

    ////fixme 混沌初始化？
    //let rand_vec_0 = numeric.random([this.dimension]);
    //let rand_vec = numeric.mul(4,rand_vec_0,numeric.sub(1,rand_vec_0));

    // 下面完成的过程是： lb + (ub - lb) * 初始值
    let range = numeric.sub(ub, lb);
    rand_vec = numeric.add(lb, numeric.mul(rand_vec, range));
    return rand_vec;
  };

  // 粒子更新公式
  this.update = function (w, c1,c2,gbest_position) {
    // 粒子更新公式为两步，分别是更新速度和位置
    // 更新速度：v(t+1) = 0.9*v(t) + 2*rand*(pbest-x) + 2*rand*(gbest-x)

    //fixme 原来为0.4-2-2或0.7298-1.494-1.494
    let velocity_part = numeric.mul(w, this.velocity);
    let pbest_part = numeric.mul(c1, numeric.random([this.dimension]), numeric.sub(this.pbest, this.position));//2
    let gbest_part = numeric.mul(c2, numeric.random([this.dimension]), numeric.sub(gbest_position, this.position));//2
    this.velocity = numeric.add(velocity_part, pbest_part, gbest_part);
    // 更新位置：x(t+1) = x(t) + v(t)
    this.position = numeric.add(this.position, this.velocity);

    //fixme 位置更新是否需要进一步约束或改进？
    // 如果更新后的位置大于ub，则放到ub上
    let gt_parts = numeric.gt(this.position, ub);
    for (let i = 0; i < gt_parts.length; i++) {
      if (gt_parts[i]) {
        //this.position[i] = this.ub[i] * (0.7 + 0.3 * Math.random());
        this.position[i] = this.ub[i];
      }
    }
    // 如果更新后的位置小于lb，则放到lb上
    let lt_parts = numeric.lt(this.position, lb);
    for (let i = 0; i < lt_parts.length; i++) {

      //混沌扰动？
      let rand0 = numeric.random([this.dimension]);
      let rand = numeric.mul(4,rand0,numeric.sub(1,rand0));

      if (lt_parts[i]) {
        //this.position[i] = this.lb[i] + 0.3 * Math.random() * (this.ub[i] - this.lb[i]);
        //this.position[i] = this.lb[i] + Math.random() * (this.ub[i] - this.lb[i]);

        this.position[i] = this.lb[i] + rand * (this.ub[i] - this.lb[i]);

      }
    }
  };

  // 求得适应度值
  this.getFitnessValue = () => {
    // 如果更新后的适应度值优于pbest值，则更新粒子历史最优
    let fitness_value = this.evaluate(this.id, this.position);
    if (fitness_value < this.pbestValue) {
      this.pbest = this.position;
      this.pbestValue = fitness_value;
    }
    return fitness_value
  }

}

export default CPSO;